home *** CD-ROM | disk | FTP | other *** search
/ Practical Algorithms for Image Analysis / Practical Algorithms for Image Analysis.iso / TARFILE.GZ / tarfile / ch_3.5 / bcd / edge / conv.c < prev    next >
Encoding:
C/C++ Source or Header  |  1999-09-11  |  4.6 KB  |  216 lines
  1. /*
  2.  * (c) Copyright 1988 by
  3.  * Robotics Principles Research Department, ATT Bell Laboratories.
  4.  * All rights reserved.
  5.  * Last modified 2/8/88 Ingemar J. Cox
  6.  * C version 8/2/88 Deborah A. Wallach
  7.  * C version 11/4/88  W.J Kropfl
  8.  */
  9. #include <stdio.h>
  10. #include <math.h>
  11. #include "edge_finder.h"
  12.  
  13. extern struct image *my_image;
  14.     static int iklim, nf, nf2;
  15.     static int iy, ny;
  16.     static int nxX2M1, ikXnx, nyX2M1Xnx, nfXnx;
  17.     static int iyMnf2Xnx, iy0, iy0Xnx, it;
  18.     static int *temp_image, *filter;
  19.     static unsigned char *tpicP;
  20.  
  21. void
  22. image_convolve(pic)
  23.      unsigned char *pic;
  24. {
  25.     register int result, ik, ix, nx;
  26.     register unsigned char *picP;
  27.     register int *filtP, *tmpP, *gaussianP;
  28.  
  29.     nx  = my_image->nx;
  30.     ny  = my_image->ny;
  31.     nf  = my_image->nf;    nfXnx = nf*nx;
  32.     nf2 = my_image->nf2;
  33.     nxX2M1 = nx + nx - 1;
  34.     nyX2M1Xnx = (ny + ny - 1)*nx;
  35.     filter = my_image->filter;
  36.  
  37. /*
  38.  *    create temp and output arrays
  39.  */
  40.  
  41.     if((temp_image = (int *)malloc(ny*nx*sizeof(int)))==0)
  42.     {
  43.         fprintf(stderr, "error: cannot allocate temp array\n");
  44.         exit(1);
  45.     }
  46.  
  47.     if((my_image->gaussian = (int *)malloc(ny*nx*sizeof(int)))==0)
  48.     {
  49.         fprintf(stderr, "error: cannot allocate gaussian array\n");
  50.         exit(1);
  51.     }
  52. /*
  53.  *    first convolve in x-dirn
  54.  *        do left side first
  55.  *        then middle
  56.  *        then right side
  57.  */
  58.     
  59.  
  60.     tmpP = &temp_image[0];
  61.     picP = pic;
  62.     for(iy=0; iy<ny; iy++, tmpP+=nx, picP+=nx)
  63.     {
  64.         for(ix=0; ix<nf2; ix++)
  65.         {
  66.             filtP = filter;
  67.             result = 0;
  68.             ik = ix - nf2;
  69.             iklim = ik + nf;
  70.             for(; ik<iklim; ik++)
  71.                 result += (*filtP++) * picP[abs(ik)];
  72.  
  73.             *(tmpP+ix) = result;
  74.         }
  75.     }
  76.     
  77.     tmpP = &temp_image[0];
  78.     tpicP = pic;
  79.     for(iy=0; iy<ny; iy++, tmpP+=nx, tpicP+=nx)
  80.     {
  81.         int xlim = nx-nf2-1;
  82.         for(ix=nf2; ix<xlim; ix++)
  83.         {
  84.             picP = tpicP+ix-nf2;
  85.             filtP = filter;
  86.             result = 0;
  87.             ik=nf;
  88.             while(ik--)
  89.                 result += (*filtP++) * (*picP++);
  90.  
  91.             *(tmpP+ix) = result;
  92.         }
  93.     }
  94.     
  95.     tmpP = &temp_image[0];
  96.     picP = pic;
  97.     for(iy=0; iy<ny; iy++, tmpP+=nx, picP+=nx)
  98.     {
  99.         for(ix=nx-nf2-1; ix<nx; ix++)
  100.         {
  101.             filtP = filter;
  102.             result = 0;
  103.             ik = ix - nf2;
  104.             iklim = ik + nf;
  105.             for(; ik<iklim; ik++)
  106.             {
  107.     /*            register int index = ix-nf2+ik;
  108.      *            if(index>=nx) index = nx+nx-index-1;
  109.      *            result += (*filtP++) * picP[index];
  110.      */
  111.                 if(ik>=nx)         /* nxX2M1 = nx+nx-ik-1; */
  112.                     result += *filtP * picP[nxX2M1-ik];
  113.                 else
  114.                     result += *filtP * picP[ik];
  115.                 filtP++;
  116.             }
  117.             *(tmpP+ix) = result;
  118.         }
  119.     }
  120.     
  121. /*
  122.  * then in y-dirn
  123.  *        first top
  124.  *        then middle
  125.  *        then bottom
  126. */
  127.     it = -nf2*nx;
  128.     for(ix=0; ix<nx; ix++)
  129.     {
  130.         gaussianP = &my_image->gaussian[ix];
  131.         tmpP = &temp_image[ix];
  132.         iy = 0; iyMnf2Xnx = it;        /* = (iy - nf2)*nx; */
  133.         for(; iy<nf2; iy++, iyMnf2Xnx+=nx)
  134.         {
  135.             filtP = filter;
  136.             result = 0;
  137.     /*        ik = iy - nf2;
  138.      *        iklim = ik + nf;
  139.      *        for(; ik<iklim; ik++)
  140.      *            result += (*filtP++) * tmpP[abs(ik)*nx];
  141.      */
  142.             ik = iyMnf2Xnx;
  143.             iklim = ik + nfXnx;
  144.              for(; ik<iklim; ik+=nx)
  145.                 result += (*filtP++) * tmpP[abs(ik)];
  146.  
  147.     /*        my_image->gaussian[iy*nx+ix] = result/my_image->norm; */
  148.             *(gaussianP) = result/my_image->norm;
  149.             gaussianP += nx;
  150.         }
  151.     }
  152.  
  153.     iy0 = nf2;
  154.     iy0Xnx = iy0 * nx;
  155.     for(ix=0; ix<nx; ix++)
  156.     {
  157.         int *col_tmpP = &temp_image[ix];
  158.         int ylim = ny-nf2-1;
  159.         iy = iy0; gaussianP = &my_image->gaussian[iy0Xnx+ix];
  160.         iyMnf2Xnx = 0;            /* = (iy-nf2)*nx; */
  161.         for( ; iy<ylim; iy++, iyMnf2Xnx+=nx)
  162.         {
  163.             filtP = filter;
  164.             tmpP = col_tmpP + iyMnf2Xnx;
  165.             result = 0;
  166.             ik=nf;
  167.             while(ik--)
  168.             {
  169.                 result += (*filtP++) * (*tmpP);
  170.                 tmpP += nx;
  171.             }
  172.     /*        my_image->gaussian[iy*nx+ix] = result/my_image->norm; */
  173.             *(gaussianP) = result/my_image->norm;
  174.             gaussianP += nx;
  175.         }
  176.     }
  177.  
  178.     iy0 = ny-nf2-1;
  179.     iy0Xnx = iy0*nx;
  180.     it = (iy0-nf2)*nx;
  181.     for(ix=0; ix<nx; ix++)
  182.     {
  183.         tmpP = &temp_image[ix];
  184.         iy = iy0; gaussianP = &my_image->gaussian[iy0Xnx+ix];
  185.         iyMnf2Xnx = it;            /* = (iy-nf2)*nx; */
  186.         for( ; iy<ny; iy++, iyMnf2Xnx+=nx)
  187.         {
  188.             filtP = filter;
  189.             result = 0;
  190.             ik = iy - nf2;
  191.             iklim = ik + nf;
  192.             ikXnx = iyMnf2Xnx;
  193.             for(; ik<iklim; ik++)
  194.             {
  195.         /*        register int index = iyMnf2+ik;
  196.          *        if(index>=ny) index = ny+ny-index-1;
  197.          *        result += (*filtP++) * tmpP[index*nx];
  198.         */
  199.                 if(ik>=ny){     /* nyX2M1Xnx = (ny+ny-1)*nx; */
  200.                     result += *filtP * tmpP[nyX2M1Xnx-ikXnx];
  201.                 }else
  202.                     result += *filtP * tmpP[ikXnx];
  203.                 filtP++;
  204.                 ikXnx += nx;
  205.             }
  206.     /*        my_image->gaussian[iy*nx+ix] = result/my_image->norm; */
  207.             *(gaussianP) = result/my_image->norm;
  208.             gaussianP += nx;
  209.         }
  210.     }
  211.     free(temp_image);
  212. #ifdef DEBUG
  213.     image_Write_int("gauss", my_image->gaussian);
  214. #endif
  215. }
  216.